NMR shifts for polycyclic aromatic hydrocarbons from first-principles 
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We present first-principles, density-functional theory calculations of the NMR chemical shifts for 
polycyclic aromatic hydrocarbons, starting with benzene and increasing sizes up to the one- and 
two-dimensional infinite limits of graphene ribbons and sheets. Our calculations are performed 
using a combination of the recently developed theory of orbital magnetization in solids, and a 
novel approach to NMR calculations where chemical shifts are obtained from the derivative of the 
orbital magnetization with respect to a microscopic, localized magnetic dipole. Using these methods 
we study on equal footing the 1 H and 13 C shifts in benzene, pyrene, coronene, in naphthalene, 
anthracene, naphthacene, and pentacene, and finally in graphene, graphite, and an infinite graphene 
ribbon. Our results show very good agreement with experiments and allow us to characterize the 
trends for the chemical shifts as a function of system size. 

PACS numbers: 71.15.-m, 71.15.Mb, 75.20.-g, 76.60.Cq 
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I. INTRODUCTION 

The experimental technique of nuclear magnetic reso- 
nance (NMR) has been known since the 1930's, 1 and has 
since become one of the most widely used methods in 
structural chemistry. 2 Along with experimental advance- 
ments over the several past decades, it was soon real- 
ized that ab-initio calculations could greatly aid in un- 
ambiguously determining structures from NMR spectra. 
Such calculations were first developed in the quantum- 
chemistry community; 3 however, these developments ap- 
plied only to finite systems, such as molecules and clus- 
ters. Application to extended systems is hindered by the 
difficulty of including macroscopic magnetic fields, which 
require a non-periodic vector potential and therefore de- 
stroy Bloch symmetry. This shortcoming was overcome 
by Mauri et ai, who developed a linear-response ap- 
proach for calculating NMR shieldings in periodic crys- 
tals based on the long-wavelength limit of a periodic 
modulation of the applied magnetic field. 4 Similarly, Se- 
bastiani and Parrinello used the Wannier representation 
to derive an alternative linear-response approach based 
on the application of an infinitesimal uniform magnetic 
field. 5 In recent years these approaches have been ad- 
vanced and refined, 6-8 leading to a growing use in modern 
plane-wave pseudopotential codes. 9-11 Nevertheless, the 
linear-response framework central to all these approaches 
makes them fairly complex and difficult to implement. 

We have recently developed a novel method to cal- 
culate chemical shifts in periodic systems, 12 where the 
need for a linear-response framework is circumvented by 
construction. While early results for simple molecules 
and hydrocarbons provided a proof of principle for our 
approach, 12 albeit limited to hydrogen shifts, we apply it 
here systematically to calculate chemical shifts for both 
hydrogen and carbon in polycyclic aromatic hydrocar- 
bons (PAH). Starting from benzene, we will explore the 
one-dimensional and two-dimensional progression that 
include either pyrene and coronene on one side, or naph- 



thalene, anthracene, etc., until we reach the infinite limits 
represented by graphene (and, by extension, graphite), or 
graphene ribbons. 

We focus this study on PAHs because of the recent 
surge of interest in energy related materials and their 
impact on environmental and health issues. PAHs occur 
in fossil fuels such as oil and coal, and they are a byprod- 
uct of incomplete combustion in e.g. wood, fat, tobacco, 
or incense. Depending on structure, these pollutants can 
be extremely toxic, carcinogenic, mutagenic, and terato- 
genic. It is thus not surprising that in recent years NMR 
techniques have been developed to determine the content 
of PAHs in our environment. 13 A tool combining ab-initio 
information with experimental techniques might lead to 
more efficient and precise devices for detecting PAHs in 
e.g. soil or air. 

This paper is organized in the following way: For 
completeness we present in Sec. II a short theory sum- 
mary of our converse approach to NMR shifts. Details 
about our numerical calculations and about the practi- 
cal implementation of the converse method can be found 
in Sec. III. Results for finite PAHs arc presented in 
Sec. IV A, whereas results for periodic systems are col- 
lected in Sec. IV B. We will conclude and summarize in 
Sec. V. 



II. THEORY 

The usual approach to calculating NMR shifts in pe- 
riodic systems is to apply a magnetic field and calcu- 
late the local field at the nucleus using a linear-response 
framework. We will refer to this approach as direct, since 
it computes the shifts directly from the applied and in- 
duced fields. Since a constant field is not compatible with 
periodic-boundary conditions, the approach developed by 
the solid-state community has been to use linear-response 
theory in the limit of long- wavelength perturbations. 4 We 
recently argued that it is actually possible to calculate the 
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shifts in periodic systems without using a linear-response 
framework. In our converse approach we circumvent the 
linear-response framework in that we relate the shifts 
to the macroscopic magnetization induced by magnetic 
point dipoles placed at the nuclear sites of interest. We 
report here the essential result and refer the reader to 
Rcf. [12] for details. 

We define E as the energy of a virtual magnetic dipolc 
m s at one nuclear center r s in the field B for a finite sys- 
tem, or as the energy per cell of a periodic lattice of such 
dipoles. Then, writing the macroscopic magnetization as 
Mp = — dE/dBp (where fi is the cell volume), 



d dE 



d dE 



dBp dm s a dm s a dBp dm s a 

'(1) 

It follows that <7 S accounts for the shielding contribution 
to the macroscopic magnetization induced by a magnetic 
point dipole m s sitting at nucleus r s and all of its pe- 
riodic replicas. Instead of applying a constant (or long- 
wavelength) field B oxt to an infinite periodic system and 
calculating the induced field at all equivalent s nuclei, we 
apply an infinite array of magnetic dipoles to all equiv- 
alent sites s, and calculate the change in magnetization. 
Since the perturbation is now periodic, it can simply be 
computed using finite differences of ground-state calcu- 
lations. Note that M = m s /0 + M ind , where the first 
term is present merely because we have included mag- 
netic dipoles by hand. The shielding is related to the true 
induced magnetization via o s , a f] = —il dM^ d / dm s ^ a . 

NMR is a technique that probes for properties of the 
electronic structure near the nuclei. Thus, a good de- 
scription of the wave functions near the nuclei is vital. 
The formalism described above can directly be imple- 
mented into all-electron programs such as LAPW codes. 
However, many codes use plane waves as a basis and the 
ionic potential is replaced by a pseudopotential to keep 
the computational cost low. In these cases, extra care 
has to be taken when calculating NMR shifts. While hy- 
drogen shifts can be described correctly using Coulombic 
(pseudo)potentials, a reconstruction of the pscudovalence 
wave functions in the core region has to be performed 
for elements beyond the first row. 4 This can be achieved 
with the so called gauge-including projector augmentcd- 
wave (GIPAW) method, which has been developed by 
Pickard and Mauri in 2001 and proved very efficient for 
NMR calculations using pseudopotentials. 7 We have de- 
rived the corresponding GIPAW formalism for the con- 
verse method, but the lengthy, mathematical details will 
be presented in a forthcoming article. 14 



III. IMPLEMENTATION AND 
CALCULATIONAL DETAILS 

We have implemented the converse method outlined 
above into the plane-wave density functional theory code 
PWSCF, which is part of the Quantum-ESPRESSO 



distribution. 15 We added an extra term to the Hamilto- 
nian taking into account the electron orbital interaction 
with a nuclear magnetic dipole m s sitting at nucleus r s , 
by the usual substitution for the momentum operator 
p — > p + | A, where A is the vector potential of a peri- 
odic array of nuclear dipoles. This is done conveniently 
in reciprocal space and requires only a few dozen lines 
of additional code. The calculation and implementation 
of the orbital magnetization 16-20 necessary to evaluate 
Eq. (1) is somewhat more involved. However, a growing 
number of codes has the calculation of the orbital mag- 
netization already implemented so that the NMR shifts 
can be calculated readily. Also, if one is only interested 
in finite systems, one can calculate the induced magnetic 
moment instead of the orbital magnetization, which is 
much simpler (the moment can be calculated from the 
quantum- mechanical current, which requires only knowl- 
edge of the wave functions). Since we are interested in 
shifts of atoms other than hydrogen, we also implemented 
a GIPAW reconstruction, as mentioned above. 

It might appear that the converse method is computa- 
tionally demanding, since we need to perform 3N cal- 
culations to obtain the shielding tensor for N atoms. 
However, note that in practice a single full self-consistent 
ground-state calculation is performed once, and then the 
dipole perturbations according to each shielding are ap- 
plied to this ground state. Convergence after the pertur- 
bation is much faster than the ground-state calculation, 
and we find that NMR shifts for systems with hundreds 
of atoms can be easily calculated. 

For our calculations in the following sections we have 
used a PBE exchange-correlation functional and norm- 
conserving pseudopotentials of the Troullier-Martins 
type 21 with a cutoff of 80 Rydberg. The structural pa- 
rameters of all systems have been optimized within this 
framework. For the dipole perturbation we chose a value 
of |m s | of 1/ib, although the results are independent of 
the exact value over a wide range. 

In order to compare the accuracy of our converse 
method to the direct method, we performed test calcu- 
lations on benzene and diamond. For the finite system 
benzene using the converse method we find an absolute 
hydrogen shielding of 22.97 ppm and carbon shielding of 

— 165.2 ppm. These results compare well with the re- 
sults obtained using the direct method: 22.69 ppm and 

— 163.4 ppm. For periodic diamond we find a shielding of 
—63.89 ppm compared to the direct method which gives 
—65.85 ppm. 6 Note that in all cases the carbon shield- 
ing includes only the shielding from the valence electrons 
and does not include the core shielding, which we calcu- 
late to be +200.3 ppm. The small deviations between 
the converse and the direct method are most likely due 
to the fact that the direct method uses a linear-response 
framework in which a finite q-vector is used to make the 
magnetic field periodic and modulate it over space, in- 
stead of using a truly constant magnetic field. 
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FIG. 1: Schematic model of benzene (left) and pyrene (right). 
In order to refer to certain atoms and their NMR shift, we use 
the labeling shown. In general, hydrogen atoms are label by 
letters, whereas carbon atoms are labeled by numerals. 




FIG. 2: Schematic model of coronene. 



IV. RESULTS 

NMR experiments usually report the isotropic shield- 
ing <j s = j;Tr[a s ] via the chemical shift S s = —((T s —a Te f). 
Here cr rc f is the isotropic shielding of a reference com- 
pound such as tetramethylsilane (TMS). Note that we 
have not calculated the shielding of TMS itself, however, 
we can still report shifts with respect to TMS by using 
our calculations for e.g. benzene as an intermediate ref- 

prpnpp- r\ — —I rr — ^benzene _ <>bcnzcnc\ p nr fbenzene 

crencc. o s — {a s er calc "tms )■ r0T °tms 

we have used the experimental values of 128.6 ppm for 
carbon and 7.26 ppm for hydrogen. 22 Henceforth, we will 
use the notation (5™ olocul ° to describe the shift relative to 
TMS of atom s in a certain molecule. We will use letters 
s = A, B, C, . . . to refer to hydrogen shifts and numerals 
s = 1,2,3,... for carbon shifts. 



A. Finite systems 

We start out by calculating the 1 H and 13 C NMR shifts 
of several small PAHs. We start with benzene, and move 
towards the one-dimensional limit by adding one ring at 
a time, to study naphthalene, anthracene, naphthacene, 
and pentacene. The two-dimensional limit is considered 
via pyrene and coronene, i.e. surrounding benzene by 
more rings in the plane. For all these calculations we 
used a large supercell with at least 16 Bohr of vacuum 
between periodic replicas. Note that the magnetic sus- 
ceptibility of these systems effectively vanishes and our 
computed shifts can be directly compared to the experi- 
mental ones without any shape correction. 23 

Our results for pyrene and coronene are collected in Ta- 
ble I. For a definition of their atomic positions see Figs. 1 
and 2. In general, we find very good agreement between 
the experimental results and our calculations. The small 
deviations are most likely due to the approximative na- 
ture of DFT with respect to the exchange and correlation 
functional. Using e.g. the local-density approximation for 
the exchange-correlation functional, we get slightly dif- 
ferent results. However, a detailed and comprehensive 
study is necessary to adequately investigate the effect 



TABLE I: Hydrogen and carbon NMR chemical shifts in ppm 
for pyrene and coronene. For a definition of atomic positions 
see Figs. 1 and 2. Note that we have used the shifts of benzene 
as a reference. 



Molecule 


atom 


experiment 


calculation 


Benzene 


1 


128.6 a 






A 


7.26 a 




Pyrene 


1 


125.5 a 


122.2 




2 


124.6 a 


122.1 




3 


130.9 a 


130.4 




4 


127.0 a 


126.0 




5 


124.6 a 


123.4 




A 


7.60-8.15 6 


8.61 




B 


7.60-8.15 b 


8.57 




C 


7.60-8.15 6 


8.47 


Coronene 


1 


126.2 a 


123.5 




2 


128.7 a 


129.3 




3 


122.6 a 


122.6 




A 


8.89 c 


9.78 



a Reference [22]. 
6 Reference [24]. 
c Reference [25]. 



of the functional used on the NMR shifts. 26 In principle, 
since we are calculating the shift through the orbital mag- 
netization, which is closely related to electrical currents, 
it might be more appropriate to use a current-density 
functional rather than one of the standard density-only 
functionals. Unfortunately, appropriate current-density 
functionals (e.g. Ref. [27]) are not yet available in most 
DFT codes. Note that the relative shifts (that is, the 
difference in shift between different atoms) are closely 
reproduced. This is important, as relative shifts play the 
key role in linking NMR data to a structure. 

As one would expect, all shifts of pyrene and coronene 
are significantly different from the pure benzene shifts. 
However, it is interesting to see that the hydrogen shifts 
increase, while all carbon shifts decrease — except one. It 
is £P yrene and S^ oroncnc that rise above the 128.6 ppm of 
benzene. These shifts both correspond to carbon atoms 
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FIG. 3: (left) Naphthalene, (right) Anthracene. 



TABLE II: Hydrogen and carbon NMR chemical shifts in ppm 
for naphthalene and anthracene. For a definition of atomic 
positions see Fig. 3. 



Molecule 


atom 


experiment 


calculation 


Benzene 


1 


128.6 a 






A 


7.26 a 




Naphthalene 


1 


125.7 c 


123.0 




2 


127.8 c 


126.2 




3 


133.4 c 


133.1 




A 


7.41 b 


8.08 




B 


7.79 6 


8.13 


Anthracene 


1 


126.1 c 


125.0 




2 


128.1 c 


127.8 




3 


131.6 c 


132.0 




4 


125.3 c 


121.6 




A 


7.41 b 


7.95 




B 


7.95 b 


8.24 




C 


8.40 6 


8.44 



a Reference [22]. 
6 Reference [24]. 
c Reference [25]. 



that have no hydrogens attached. 

Our results for the one-dimensional chains from naph- 
thalene to pentacene are given in Tables II and III. For 
a definition of the corresponding atomic positions see 
Figs. 3 and 4. As before, we find very good agreement be- 
tween the experimental results and our calculations. The 
exact same statements concerning the differences com- 
pared to benzene apply: hydrogen shifts increase, while 
carbon shifts decrease — expect one. In these chains it is 
the carbon shift of atom 3 that increases, which again 
corresponds to a carbon that is not connected to a hy- 
drogen. 

Note that we have not included experimental results 
for naphthacene and pentacene. Although there are some 
results reported (see e.g. Ref. [28]) these molecules are 
insoluble in most of the solvents used for the preparation 
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FIG. 4: (left) Naphthacene. (right) Pentacene. 

TABLE III: Hydrogen and carbon NMR chemical shifts in 
ppm for naphthacene and pentacene. For a definition of 
atomic positions see Fig. 4. 

Naphthacene Pentacene 



atom 


calc. 


atom 


calc. 


1 


122.3 


1 


123.3 


2 


126.5 


2 


127.4 


3 


130.5 


3 


130.6 


4 


125.5 


4 


125.0 


5 


128.1 


5 


128.5 


A 


8.31 


6 


125.4 


B 


8.68 


A 


8.37 


C 


9.07 


B 


8.74 






C 


9.21 






D 


9.46 



of NMR cuvettes. 



B. Periodic systems 

The systems discussed above can systematically be ex- 
panded until they eventually become infinite. In a one- 
dimensional way, we extend the chains naphthalene, an- 
thracene, naphthacene, pentacene, . . . until we reach an 
infinite ribbon. On the other hand, if we extend pyrene 
and coronene in a two-dimensional way, we arrive at 
graphene. Stacking graphene sheets on top of each other 
yields graphite. 

We have calculated the hydrogen and carbon shifts 
for these extended systems. These infinite systems show 
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FIG. 5: Section of an infinite ribbon consisting of benzene 
rings. 

metallic behavior and therefore require a dense k-point 
mesh. For our calculations we have used up to 16 k- 
points in the directions of periodicity and a gaussian 
smearing of 0.001 Rydberg to obtain converged results. 
As another result of the metallic behavior, Knight shifts 
might influence the NMR spectrum. While schemes do 
exist to calculate Knight shifts within plane-wave DFT 
approaches, 29 we did not include them in our calcula- 
tions. 

For graphene we find a carbon shift of 118.0 ppm. 
The two in-equivalent carbon shifts in graphite give 
124.3 ppm and 134.9 ppm. For the infinite rib- 
bon depicted in Fig. 5 we find (5^ bbon =8.56 ppm, 
£ribbon =128 o ppnlj and <jribbon =132 2 ppm . Experimen- 
tally, solid state NMR spectra of graphite show a broad 
peak in the range 155-^179 ppm, depending on sample 
preparation. This broadening is due to a fairly long lat- 
tice relaxation time T\ and the presence of conduction 
electrons. 30 

In order to compare our results for graphite with ex- 
perimental data, in principle one would have to include a 
shape correction involving the susceptibility. 23 This cor- 
rection is estimated to be of the order of 1-^2 ppm in 
most forms of graphite sample such as powder samples. 
The only exception is highly oriented pyrolytic graphite, 
which has a very large and anisotropic magnetic suscep- 
tibility, where the shape correction is expected to be of 
the order of 10-f-20 ppm. 

The calculated shifts for these periodic systems are in- 



teresting in and of themselves. However, we now want to 
focus on the change in shift as the systems grow larger 
and eventually become infinite. One would expect to 
find a correlation between the shifts in finite and infi- 
nite systems. Indeed, looking at the carbon shifts we 

find <5 bcnzcnc — > ^Pyrcne _^ ^oroncne Jgraphcnc w j th 

the results 128.6 -> 123.4 -> 122.6 -> 118.0 ppm. Al- 
though NMR is a local probe, it turns out that the size 
of coronene is not yet big enough for its bulk atoms to 
behave like graphene. 

If we start at benzene again and add rings in a liner 
fashion, we can grow a chain of arbitrary length. Looking 
at the carbon shifts we find the sequence ^| nthracene _ ► 
^pcntaccnc _^ gibbon with thc corresponding results 
121.6 — > 125.4 — > 128.0 ppm. As with the graphene 
sheet, this result suggests that a finite chain would need 
to be longer than pentacene so that its carbon atoms be- 
have like an infinite ribbon. Looking at sequences that 
end at <^ bbon or #j ibbon , no clear trend is visible. 

V. CONCLUSIONS 

We have used an alternative first-principles method to 
calculate NMR chemical shifts of a variety of polycyclic 
aromatic hydrocarbons and related infinite systems. Our 
results are in good agreement with experiment. By going 
from finite to periodic systems, we observe trends in shifts 
which suggest that neither coronene nor pentacene are 
large enough to model their infinite counterparts. 

In future work we would like to include thc shifts 
of other PAHs as well as study the Knight shifts of 
graphene. 
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